Fase 3. Modelación

Exportemos a R nuestro conjunto de datos para poder comenzar el análisis:

# Exportamos los datos con la función import de la 
# librería rio.

library(readxl)
## Warning: package 'readxl' was built under R version 4.3.2
datos_sucios<- read_excel("AirQualityUCI.xlsx")

Para mejorar la claridad y evitar confusiones, vamos a proceder a renombrar todas las variables.

# Renombramos las variables para que sea más fácil poder 
# trabajar con ellas.

library(dplyr)
## Warning: package 'dplyr' was built under R version 4.3.2
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
datos_sucios<- rename(datos_sucios, 
                Fecha = "Date",
                Hora = "Time",
                Monoxido_carbono = "CO(GT)",
                Oxido_estaño = "PT08.S1(CO)",
                Hidrocarburos_no_metanicos = "NMHC(GT)",
                Benceno = "C6H6(GT)",
                Titanio = "PT08.S2(NMHC)",
                Oxidos_nitrogeno = "NOx(GT)",
                Oxido_tungsteno = "PT08.S3(NOx)",
                Dioxido_nitrogeno = "NO2(GT)",
                Oxido_tungsteno_NO2 = "PT08.S4(NO2)",
                Oxido_indio = "PT08.S5(O3)",
                Temperatura = "T",
                Humedad_relativa  = "RH",
                Humedad_absoluta = "AH")

Comprobemos si el cambio de nombres se realizó de forma adecuada:

# Veamos si se realizó el cambio en los nombres: 
names(datos_sucios)
##  [1] "Fecha"                      "Hora"                      
##  [3] "Monoxido_carbono"           "Oxido_estaño"              
##  [5] "Hidrocarburos_no_metanicos" "Benceno"                   
##  [7] "Titanio"                    "Oxidos_nitrogeno"          
##  [9] "Oxido_tungsteno"            "Dioxido_nitrogeno"         
## [11] "Oxido_tungsteno_NO2"        "Oxido_indio"               
## [13] "Temperatura"                "Humedad_relativa"          
## [15] "Humedad_absoluta"

Con las variables renombradas y una nomenclatura más clara, veamos un resumen rápido de los datos:

# Resumen rapido, conociendo los tipos de datos 
skimr::skim(datos_sucios)
Data summary
Name datos_sucios
Number of rows 9357
Number of columns 15
_______________________
Column type frequency:
numeric 13
POSIXct 2
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Monoxido_carbono 0 1 -34.21 77.66 -200 0.60 1.50 2.60 11.90 ▂▁▁▁▇
Oxido_estaño 0 1 1048.87 329.82 -200 921.00 1052.50 1221.25 2039.75 ▁▁▇▅▁
Hidrocarburos_no_metanicos 0 1 -159.09 139.79 -200 -200.00 -200.00 -200.00 1189.00 ▇▁▁▁▁
Benceno 0 1 1.87 41.38 -200 4.00 7.89 13.64 63.74 ▁▁▁▇▅
Titanio 0 1 894.48 342.32 -200 711.00 894.50 1104.75 2214.00 ▁▅▇▂▁
Oxidos_nitrogeno 0 1 168.60 257.42 -200 50.00 141.00 284.20 1479.00 ▇▆▂▁▁
Oxido_tungsteno 0 1 794.87 321.98 -200 637.00 794.25 960.25 2682.75 ▁▇▃▁▁
Dioxido_nitrogeno 0 1 58.14 126.93 -200 53.00 96.00 133.00 339.70 ▃▁▇▅▁
Oxido_tungsteno_NO2 0 1 1391.36 467.19 -200 1184.75 1445.50 1662.00 2775.00 ▁▂▇▅▁
Oxido_indio 0 1 974.95 456.92 -200 699.75 942.00 1255.25 2522.75 ▁▇▇▃▁
Temperatura 0 1 9.78 43.20 -200 10.95 17.20 24.07 44.60 ▁▁▁▁▇
Humedad_relativa 0 1 39.48 51.22 -200 34.05 48.55 61.88 88.73 ▁▁▁▂▇
Humedad_absoluta 0 1 -6.84 38.98 -200 0.69 0.98 1.30 2.23 ▁▁▁▁▇

Variable type: POSIXct

skim_variable n_missing complete_rate min max median n_unique
Fecha 0 1 2004-03-10 2005-04-04 00:00:00 2004-09-21 00:00:00 391
Hora 0 1 1899-12-31 1899-12-31 23:00:00 1899-12-31 11:00:00 24

Con el resumen podemos observar que el conjunto de datos cuenta con 9357 observaciones de 15 variables. R ha identificado dos tipos de datos en el conjunto: numeric (15 variables) y POSIXct (2 variables). Las variables de tipo numeric corresponden a las mediciones de los sensores, mientras que las dos variables de tipo POSIXct son la fecha y hora en que se realizaron las mediciones.
Aunque parece que los datos están completos, es importante recordar que los valores faltantes en el conjunto fueron etiquetados con un valor de -200. Por lo tanto, es necesario determinar la cantidad real de datos faltantes en el conjunto tomando en cuenta esta etiqueta.

# Veamos un conteo de los datos faltantes para cada
# una de las variables.

sapply(datos_sucios, function(x) sum(x == -200))
##                      Fecha                       Hora 
##                          0                          0 
##           Monoxido_carbono               Oxido_estaño 
##                       1683                        366 
## Hidrocarburos_no_metanicos                    Benceno 
##                       8443                        366 
##                    Titanio           Oxidos_nitrogeno 
##                        366                       1639 
##            Oxido_tungsteno          Dioxido_nitrogeno 
##                        366                       1642 
##        Oxido_tungsteno_NO2                Oxido_indio 
##                        366                        366 
##                Temperatura           Humedad_relativa 
##                        366                        366 
##           Humedad_absoluta 
##                        366

Ahora que hemos evaluado la cantidad real de datos faltantes en el conjunto de datos, podemos observar que la variable que contiene los datos de la concentración promedio de hidrocarburos no metánicos es la que tiene el mayor número de datos faltantes, con un total de 8443. Además, hay tres variables con más de 1600 datos faltantes y finalmente 9 variables con 366 datos faltantes, que corresponden a los datos obtenidos del sensor.

Para poder utilizar algunas de las funciones útiles de R destinadas específicamente a tratar con datos faltantes, es necesario cambiar los valores de -200 por NA. De esta manera, podremos aprovechar al máximo las herramientas disponibles en R para el análisis de datos.

# Reemplzamos los -200 por NA
datos_sucios <- replace(datos_sucios, datos_sucios == -200, NA)

# Usamos la libreria naniar para ver de forma grafica
# información de los datos faltantes
library(naniar)
## Warning: package 'naniar' was built under R version 4.3.2
# Porcentaje de NA
gg_miss_var(datos_sucios, show_pct = T)

# ¿En dónde se localizan los datos faltantes?:
vis_miss(datos_sucios)

De todas las variables disponibles en nuestro conjunto de datos nos centraremos en las siguientes variables, las cuales son fundamentales para la evaluación de la calidad del aire y son ampliamente monitoreadas en diferentes lugares del mundo:

  1. Concentración promedio de \(CO\).
  2. Concentración promedio de benceno.
  3. Concentración promedio de \(NO_x\).
  4. Concentración promedio de \(NO_2\).
  5. Temperatura.
  6. Humedad relativa.

Aunque las concentraciones de los compuestos químicos en la atmósfera son lo más relevante, consideramos la temperatura y la humedad relativa importantes para la evaluación de la calidad del aire, ya que pueden afectar la formación y dispersión de contaminantes en el aire.

Después de ver los datos faltantes y ver que la variable de la concentración promedio de hidrocarburos no metánicos es la que más datos faltantes tiene, la haremos a un lado para realizar la modelación con la mayor cantidad de datos.

# Seleccionamos las variables que usaremos para el 
# análisis

library(dplyr)
calidad_Aire<- select(datos_sucios, 
                      Monoxido_carbono, 
                      Benceno,
                      Oxidos_nitrogeno,
                      Dioxido_nitrogeno,
                      Temperatura,
                      Humedad_relativa)

Finalmente, eliminemos todos los renglones donde exista al menos un dato faltante de nuestro conjunto de datos.

# Eliminamos renglones con CUALQUIER valor faltante

library(tidyr)
## Warning: package 'tidyr' was built under R version 4.3.2
calidad_Aire_Limpio<- drop_na(calidad_Aire)

paste("La cantidad de observaciones que quedan después de eliminar los datos faltantes son: ", nrow(calidad_Aire_Limpio))
## [1] "La cantidad de observaciones que quedan después de eliminar los datos faltantes son:  6941"

Vemos que después de eliminar los renglones con al menos un dato faltante nos quedan nuestras 6 variables con 6941 observaciones; es decir, eliminamos aproximadamente el 26% de los datos.

Consideremos las variables que hemos seleccionado, es importante estudiar sus distribuciones de probabilidad. Recordando un poco lo que hicimos en la Fase 2 veamos si estas cumplen o no con tener un perfil “normal”, veamos las distribuciones:
NOTA: La línea roja es la curva normal teórica y la línea azul es la densidad empírica.

library(psych)
## Warning: package 'psych' was built under R version 4.3.2
library(paletteer) # Librería para las paletas de colores
## Warning: package 'paletteer' was built under R version 4.3.2
# Histogramas

multi.hist(x = calidad_Aire_Limpio, bcol=paletteer_d("RColorBrewer::BuPu"), dcol = c("blue", "red"), dlty = c("dotted", "solid"),lwd=c(2,3),
           main = c("Monóxido de Carbono","Benceno",
                    "Óxidos de Nitrógeno","Dióxido de Nitrógeno",
                    "Temperatura","Humedad Relativa"), global = FALSE)

A primera impresión notamos que las variables correspondientes a las concentraciones de monóxido de carbono, benceno y óxidos de nitrógeno no parecen tener este perfil normal que estamos buscando, mientras que las variables que corresponden a la concentración de dióxido de nitrógeno, temperatura y humedad relativa parece que si cumplen con el perfil normal. Pero veamos unas graficas más:

attach(calidad_Aire_Limpio)

par(mfrow=c(2,3))

qqnorm(Monoxido_carbono, main="Normal Q-Q Plot (Monóxido de Carbono)",col="darkblue")
qqline(Monoxido_carbono, col="red", lwd=2)

qqnorm(Benceno, main="Normal Q-Q Plot (Benceno)",col="darkblue")
qqline(Benceno, col="red", lwd=2)

qqnorm(Oxidos_nitrogeno, main="Normal Q-Q Plot (Óxidos de Nitrógeno)",col="darkblue")
qqline(Oxidos_nitrogeno, col="red", lwd=2)

qqnorm(Dioxido_nitrogeno, main="Normal Q-Q Plot (Dióxido de Nitrógeno)",col="darkblue")
qqline(Dioxido_nitrogeno, col="red", lwd=2)

qqnorm(Temperatura, main="Normal Q-Q Plot (Temperatura)",col="darkblue")
qqline(Temperatura, col="red", lwd=2)

qqnorm(Humedad_relativa, main="Normal Q-Q Plot (Humedad Relativa)",col="darkblue")
qqline(Humedad_relativa, col="red", lwd=2)

Una forma de ver si existe normalidad consiste en comparar los cuantiles de la distribución observada con los cuantiles teóricos de una distribución normal con la misma media y desviación estándar que los datos. Cuanto más se aproximen los datos a una normal, más alineados están los puntos entorno a la recta. Esta comparación corresponde a las graficas anteriores, con eso podemos concluir que nuestras sospechas eran ciertas, las observaciones de las variables correspondientes a las concentraciones de monóxido de carbono, benceno y óxidos de nitrógeno no caen sobre la recta lo que nos indica que no hay normalidad, para las otras tres variables podríamos decir que la “gran mayoría” de observaciones si caen sobre la recta. Incluso podríamos realizar una prueba mas rigurosa para comprobar la normalidad, pero por esta ocasión nos quedaremos y aceptaremos estas conclusiones. Entonces, tenemos tres variables problemáticas:

Profundicemos en detalle en estas tres variables problemáticas y examinemos qué distribuciones de probabilidad podrían modelarlas de manera más precisa.

A primera vista, las distribuciones de probabilidad de las tres variables problemáticas parecen exhibir un comportamiento similar (es curioso que tres variables sean descritas por una distribución normal y que las otras tres variables no tengan este comportamiento normal pero que sí sigan un comportamiento bastante parecido entre ellas). Si nos aventuramos aún más, podríamos afirmar que podemos modelar su comportamiento utilizando distribuciones Gamma. Sin embargo, sería prudente no aceptar esto como un hecho y analizar más detenidamente las variables:

Veamos primero la gráfica de Cullen & Frey, esta gráfica se representa el coeficiente de asimetría en el eje x y el coeficiente de curtosis en el eje y. Cada punto en el gráfico representa una distribución de datos específica. Si los datos siguen una distribución normal, los puntos se agruparán cerca del (0,3) en el gráfico. Si los datos tienen una asimetría o curtosis significativa, los puntos se alejarán de ese punto y se desplazarán hacia diferentes direcciones.
La grafica de Cullen & Frey nos ayuda a escoger alguno de los modelos de probabilidad basados en la curtosis y la simetría.

# Creamos el grafico Cullen & Frey
library(fitdistrplus)
## Warning: package 'fitdistrplus' was built under R version 4.3.2
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## Loading required package: survival
descdist(data = Monoxido_carbono, graph = TRUE)

## summary statistics
## ------
## min:  0.1   max:  11.9 
## median:  1.9 
## mean:  2.182467 
## estimated sd:  1.441158 
## estimated skewness:  1.33915 
## estimated kurtosis:  5.620881

Basado en el perfil del histograma y la gráfica de Cullen & Frey, nuestras observaciones preliminares parecen confirmarse, y es plausible que el modelo de probabilidad más apropiado sea una Distribución Gamma, caracterizada por los parámetros alfa y beta.
Veamos ahora el valor del AIC (Akaike’s Information Criterion) para diferentes distribuciones. El AIC es una medida utilizada para la selección de modelos en el contexto del análisis de datos y el ajuste de modelos estadísticos. Cuanto más bajo sea el valor del AIC, mejor se considera que es el ajuste del modelo. Esto significa que un modelo con un AIC más bajo tiene un buen equilibrio entre ajuste y complejidad, y se prefiere sobre modelos con AIC más alto.

library(univariateML)
## Warning: package 'univariateML' was built under R version 4.3.2
library(dplyr)
library(tibble)
## Warning: package 'tibble' was built under R version 4.3.2
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 8):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
# Calculemos el AIC para diferentes distribuciones de probabilidad
# y así poder elegir la mejor distribución.

comparacion_aic <- AIC(
  mlbetapr(Monoxido_carbono), # Beta prime distribution
  mlexp(Monoxido_carbono), # Exponential distribution
  mlinvgamma(Monoxido_carbono), # Inverse Gamma distribution
  mlgamma(Monoxido_carbono), # Gamma distribution
  mllnorm(Monoxido_carbono), # Log-normal distribution
  mlrayleigh(Monoxido_carbono), # Rayleigh distribution
  mlinvgauss(Monoxido_carbono), # Inverse Gaussian
  mlweibull(Monoxido_carbono), # Weibull distribution
  mlinvweibull(Monoxido_carbono) # Inverse Weibull distribution
)
comparacion_aic %>% rownames_to_column(var = "distribucion") %>% 
  arrange(AIC)
##                     distribucion df      AIC
## 1      mlgamma(Monoxido_carbono)  2 22420.80
## 2    mlweibull(Monoxido_carbono)  2 22602.19
## 3      mllnorm(Monoxido_carbono)  2 22810.90
## 4     mlbetapr(Monoxido_carbono)  2 22990.51
## 5   mlrayleigh(Monoxido_carbono)  1 23279.82
## 6   mlinvgauss(Monoxido_carbono)  2 23348.53
## 7   mlinvgamma(Monoxido_carbono)  2 24604.22
## 8        mlexp(Monoxido_carbono)  1 24718.29
## 9 mlinvweibull(Monoxido_carbono)  2 25344.70

Al analizar los valores del AIC de las diferentes distribuciones, se observa que el valor más bajo corresponde a la distribución Gamma, lo que indica que efectivamente es la distribución más adecuada.
Hagamos el ajuste con una distribución Gamma:

library(latex2exp)
## Warning: package 'latex2exp' was built under R version 4.3.2
library(MASS)
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 4):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 10):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
library(paletteer)

#Realizamos el ajuste de la función de densidad Gamma
monoxido_fitGamma<- fitdistr(Monoxido_carbono, densfun="gamma")

# Graficamos el histograma junto con la curva de ajuste
hist(Monoxido_carbono,col=paletteer_d("ggsci::indigo_material"),
     main=TeX("Monóxido de Carbono"),
     xlab= TeX("$mg/m^3$"),
     ylim = c(0, 0.4),
     breaks=20,
     prob=TRUE)
curve(dgamma(x, monoxido_fitGamma$estimate[1],
                monoxido_fitGamma$estimate[2]),col="red",
                lwd=2,add=T)

paste("Los parámetros encontrados con R por el método de máxima verosimilitud son: ", monoxido_fitGamma$estimate[1], monoxido_fitGamma$estimate[2])
## [1] "Los parámetros encontrados con R por el método de máxima verosimilitud son:  2.34925510077524 1.07642238920912"
Entonces, nuestra variable de concentración de Monóxido de Carbono tiene el comportamiento:

\(CO \sim Gamma(\alpha = 2.349, \beta= 1.076)\)


Para las otras dos variables problemáticas seguiremos el mismo procedimiento.

Veamos primero la gráfica de Cullen & Frey para la variable que contiene los datos de la concentración de Benceno:

# Creamos el grafico Cullen & Frey
library(fitdistrplus)
descdist(data = Benceno, graph = TRUE)

## summary statistics
## ------
## min:  0.1815254   max:  63.74148 
## median:  8.788282 
## mean:  10.55441 
## estimated sd:  7.46517 
## estimated skewness:  1.301695 
## estimated kurtosis:  5.409269

Vemos que de acuerdo al perfil del histograma y la gráfica de Cullen & Frey parece que nuestras sospechas se hacen realidad, el modelo de probabilidad que parece más adecuado parece ser nuevamente una Distribución Gamma de parámetros alfa y beta.
Veamos ahora el valor del AIC (Akaike’s Information Criterion) para diferentes distribuciones.

library(univariateML)
library(dplyr)
library(tibble)
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 3):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 5):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 11):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
# Calculemos el AIC para diferentes distribuciones de probabilidad
# y así poder elegir la mejor distribución.

comparacion_aic <- AIC(
  mlbetapr(Benceno), # Beta prime distribution
  mlexp(Benceno), # Exponential distribution
  mlinvgamma(Benceno), # Inverse Gamma distribution
  mlgamma(Benceno), # Gamma distribution
  mllnorm(Benceno), # Log-normal distribution
  mlrayleigh(Benceno), # Rayleigh distribution
  mlinvgauss(Benceno), # Inverse Gaussian
  mlweibull(Benceno), # Weibull distribution
  mlinvweibull(Benceno) # Inverse Weibull distribution
)
comparacion_aic %>% rownames_to_column(var = "distribucion") %>% 
  arrange(AIC)
##            distribucion df      AIC
## 1      mlgamma(Benceno)  2 45042.35
## 2    mlweibull(Benceno)  2 45125.36
## 3      mllnorm(Benceno)  2 45583.01
## 4   mlinvgauss(Benceno)  2 46333.64
## 5   mlrayleigh(Benceno)  1 46411.88
## 6        mlexp(Benceno)  1 46597.54
## 7     mlbetapr(Benceno)  2 46711.92
## 8   mlinvgamma(Benceno)  2 47774.08
## 9 mlinvweibull(Benceno)  2 48181.68

Al analizar los valores del AIC de las diferentes distribuciones, se observa que el valor más bajo corresponde a la distribución Gamma, lo que indica que efectivamente es la distribución más adecuada.
Hagamos el ajuste con una distribución Gamma:

library(latex2exp)
library(MASS)
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 3):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 4):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 6):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 12):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
library(paletteer)

#Realizamos el ajuste de la función de densidad Gamma
benceno_fitGamma<- fitdistr(Benceno, densfun="gamma")
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced

## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
# Graficamos el histograma junto con la curva de ajuste
hist(Benceno,col=paletteer_d("ggsci::deep_orange_material"),
     main=TeX("Benceno"),
     xlab= TeX("$\\mu g/m^3$"),
     ylim = c(0, 0.08),
     breaks=20,
     prob=TRUE)
curve(dgamma(x, benceno_fitGamma$estimate[1],
                benceno_fitGamma$estimate[2]),col="red",
                lwd=2,add=T)

paste("Los parámetros encontrados con R por el método de máxima verosimilitud son: ", benceno_fitGamma$estimate[1], benceno_fitGamma$estimate[2])
## [1] "Los parámetros encontrados con R por el método de máxima verosimilitud son:  1.97409183844 0.187039581597446"
Entonces, nuestra variable de concentración de Benceno tiene el comportamiento:

\(Benceno \sim Gamma(\alpha = 1.974, \beta= 0.187)\)


Estudiemos nuestra última variable problemática. Primero la gráfica de Cullen & Frey:

# Creamos el grafico Cullen & Frey
library(fitdistrplus)
descdist(data = Oxidos_nitrogeno, graph = TRUE)

## summary statistics
## ------
## min:  2   max:  1479 
## median:  186 
## mean:  250.6565 
## estimated sd:  208.604 
## estimated skewness:  1.643631 
## estimated kurtosis:  6.155228

Para esta última variable, parece que la Distribución Gamma ya no es la más adecuada para ajustar los datos. Parece ser que una distribución beta podría ser la mejor opción. Sin embargo, antes de sacar conclusiones, es importante evaluar el valor del AIC (Akaike’s Information Criterion) para diferentes distribuciones.

library(univariateML)
library(dplyr)
library(tibble)
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 3):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 4):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 5):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 7):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 13):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
# Calculemos el AIC para diferentes distribuciones de probabilidad
# y así poder elegir la mejor distribución.

comparacion_aic <- AIC(
  mlbetapr(Oxidos_nitrogeno), # Beta prime distribution
  mlexp(Oxidos_nitrogeno), # Exponential distribution
  mlinvgamma(Oxidos_nitrogeno), # Inverse Gamma distribution
  mlgamma(Oxidos_nitrogeno), # Gamma distribution
  mllnorm(Oxidos_nitrogeno), # Log-normal distribution
  mlrayleigh(Oxidos_nitrogeno), # Rayleigh distribution
  mlinvgauss(Oxidos_nitrogeno), # Inverse Gaussian
  mlweibull(Oxidos_nitrogeno), # Weibull distribution
  mlinvweibull(Oxidos_nitrogeno) # Inverse Weibull distribution
)
comparacion_aic %>% rownames_to_column(var = "distribucion") %>% 
  arrange(AIC)
##                     distribucion df      AIC
## 1      mllnorm(Oxidos_nitrogeno)  2 89673.24
## 2      mlgamma(Oxidos_nitrogeno)  2 89693.83
## 3    mlweibull(Oxidos_nitrogeno)  2 89884.58
## 4   mlinvgauss(Oxidos_nitrogeno)  2 90044.34
## 5        mlexp(Oxidos_nitrogeno)  1 90569.33
## 6     mlbetapr(Oxidos_nitrogeno)  2 91063.12
## 7   mlinvgamma(Oxidos_nitrogeno)  2 91137.15
## 8 mlinvweibull(Oxidos_nitrogeno)  2 91448.87
## 9   mlrayleigh(Oxidos_nitrogeno)  1 92914.62

Al analizar los valores del AIC de las diferentes distribuciones, se observa que el valor más bajo corresponde a la distribución Log-Normal, una distribución que no teníamos en mente.
Hagamos el ajuste con una distribución Log-Normal:

library(latex2exp)
library(MASS)
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 3):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 4):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 5):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 6):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 8):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 14):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
library(paletteer)

#Realizamos el ajuste de la función de densidad Gamma
oxidos_fitLogNormal<- fitdistr(Oxidos_nitrogeno, densfun="lognormal")

# Graficamos el histograma junto con la curva de ajuste
hist(Oxidos_nitrogeno,col=paletteer_d("ggsci::teal_material"),
     main=TeX("Óxidos de Nitrógeno"),
     xlab= TeX("$ppb$"),
     ylim = c(0, 0.0045),
     breaks=20,
     prob=TRUE)
curve(dlnorm(x, oxidos_fitLogNormal$estimate[1],
                oxidos_fitLogNormal$estimate[2]),col="red",
                lwd=2,add=T)

paste("Los parámetros encontrados con R por el método de máxima verosimilitud son: ", oxidos_fitLogNormal$estimate[1], oxidos_fitLogNormal$estimate[2])
## [1] "Los parámetros encontrados con R por el método de máxima verosimilitud son:  5.18820287112231 0.862644837560372"

Entonces, podemos moodelar el comportamiento de la variable de concentración de Óxidos de Nitrógeno con una distribución Log-Normal con parametros:
meanlog = 5.188, sdlog=0.862

Método de Box-Cox

Para transformar la(s) variable(s) para que tenga una distribución normal usaremos la transformación de Box-Cox que es una técnica estadística utilizada para estabilizar la varianza y mejorar la normalidad de los datos. Esta técnica fue propuesta por los estadísticos George Box y David Cox.
Para cada una de las variables problemáticas realizaremos estas transformaciones y mostraremos a continuación los histogramas y las gráficas Normal Q-Q para ver si la transformación funcionó y las variables se normalizaron.

# Mónoxido de Carbono
library(forecast)
## Warning: package 'forecast' was built under R version 4.3.2
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 4):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 5):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 6):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 7):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 8):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 10):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 16):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
library(paletteer)
library(latex2exp)

# Encontramos lambdas optimas
lambda_monoxido = BoxCox.lambda(Monoxido_carbono,method="loglik",
                                lower=-100,upper=100)
lambda_benceno = BoxCox.lambda(Benceno,method="loglik",
                               lower=-100,upper=100)
lambda_oxidos = BoxCox.lambda(Oxidos_nitrogeno,method="guerrero",
                              lower=-100,upper=100)
## Warning in optimize(guer.cv, c(lower, upper), x = x, nonseasonal.length =
## nonseasonal.length): NA/Inf replaced by maximum positive value
#Realizamos las transformaciones
trans.monoxido<- BoxCox(Monoxido_carbono, lambda_monoxido)
trans.benceno<- BoxCox(Benceno, lambda_benceno)
trans.oxidos<- BoxCox(Oxidos_nitrogeno, lambda_oxidos)

# Histogramas
# Monóxido de Carbono
par(mfrow=c(3,3))
hist(Monoxido_carbono,col=paletteer_d("ggsci::indigo_material"),
     main=TeX("Monóxido de Carbono"),
     breaks=20,
     ylim = c(0, 0.4),
     prob=TRUE) # variable original

hist(trans.monoxido, col=paletteer_dynamic("cartography::orange.pal", 20),
     main=TeX("Transformación"),
     breaks=20,
     ylim = c(0, 0.5),
     prob=TRUE)  # variable transformada
# Graficas Q-Q
qqnorm(trans.monoxido, col="darkblue", cex = 0.7)
qqline(trans.monoxido, col="red", lwd=1)

# Benceno
hist(Benceno,col=paletteer_d("ggsci::deep_orange_material"),
     main=TeX("Benceno"),
     breaks=20,
     ylim = c(0, 0.08),
     prob=TRUE) # variable original

hist(trans.benceno, col=paletteer_dynamic("cartography::orange.pal", 20),
     main=TeX("Transformación"),
     breaks=20,
     ylim = c(0, 0.3),
     prob=TRUE)  # variable transformada
# Graficas Q-Q
qqnorm(trans.benceno,col="darkblue", cex = 0.7)
qqline(trans.benceno, col="red", lwd=1)

# Óxidos de Nitrógeno
hist(Oxidos_nitrogeno,col=paletteer_d("ggsci::teal_material"),
     main=TeX("Óxidos de Nitrógeno"),
     breaks=20,
     ylim = c(0, 0.0045),
     prob=TRUE) # variable original

hist(trans.oxidos, col=paletteer_dynamic("cartography::orange.pal", 20),
     main=TeX("Transformación"),
     breaks=20,
     ylim = c(0, 0.3),
     prob=TRUE)  # variable transformada
# Graficas Q-Q
qqnorm(trans.oxidos,col="darkblue", cex = 0.7)
qqline(trans.oxidos, col="red", lwd=1)

Podemos ver que después de la transformación realizada a las tres variables problemáticas, los perfiles de sus distribuciones (transformadas) y las gráficas Q-Q nos muestran aparente normalidad en todos los casos.

El monóxido de carbono (\(CO\)) es un gas resultante de la combustión incompleta, que ocurre cuando existe una baja concentración de oxígeno. Según la bibliografía, aproximadamente el 86% de las emisiones de \(CO\) provienen del sector del transporte, seguido por un 6% atribuido a la quema de combustibles en la industria, y un 3% originado en procesos industriales. El 4% restante se genera a través de quemas y otros procesos no identificados. Asimismo, el CO se produce de manera natural mediante la oxidación del metano, un proceso comúnmente asociado con la descomposición de materia orgánica.

El \(CO\) puede tener efectos adversos en la salud, ya que compite con el oxígeno en la corriente sanguínea, disminuyendo así la capacidad de transporte de oxígeno a los diferentes órganos. Las personas sensibles, especialmente aquellas que padecen problemas cardíacos, pueden experimentar una reducción en su capacidad de oxigenación.

Dado lo anteriormente expuesto, resulta crucial llevar a cabo un monitoreo de los niveles de CO. Por tanto, hemos decidido seleccionar la concentración de monóxido de carbono (\(CO\)) como nuestra variable respuesta.

Para realizar la división del conjunto de datos en dos grupos, se propuso inicialmente basarla en la presencia de registros de concentración de monóxido de carbono que excedieran los límites establecidos por las organizaciones pertinentes. Por ejemplo, en México, el límite para la concentración de monóxido de carbono como contaminante atmosférico es de 12.595 miligramos por metro cúbico en un promedio móvil de ocho horas al año, con el objetivo de proteger la salud de la población susceptible. De manera similar, la Agencia de Protección Ambiental de Estados Unidos (EPA, por sus siglas en inglés) ha establecido un límite ambiental de 10 miligramos por metro cúbico para el monóxido de carbono en el aire, promediado en un período de ocho horas.

Sin embargo, surge un problema con esta división, ya que hay muy pocos registros que superen estas concentraciones. Si consideramos únicamente las observaciones que cumplen con un nivel superior a los 10 miligramos por metro cúbico, nos encontramos con tan solo 5 observaciones. Por lo tanto, se ha tomado la decisión de separar el conjunto de datos basándonos en la temperatura. Dado que las observaciones se obtuvieron en una ciudad italiana (no se especifica cuál), se han revisado los registros de temperatura promedio en la ciudad de Roma, Italia, correspondientes al año 2004, y se ha encontrado una temperatura media anual de 17.1°C. Al analizar los datos, se observa que el promedio de temperatura es de 17.7°C, lo que parece ser un criterio adecuado para la división. Además, resultará interesante comparar el promedio de la variable respuesta entre ambos grupos para ver si la concentración de monóxido de carbono aumenta o disminuye en función de la presencia de temperaturas por encima del promedio.

  1. Conjunto 1. Temperatura \(\leq\) 17.7°C
  2. Conjunto 2. Temperatura \(>\) 17.7°C

Hagamos entonces la separación:

attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 3):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 5):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 6):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 7):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 8):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 9):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 11):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 17):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
datos_baja_temperatura<- calidad_Aire_Limpio[Temperatura <= 17.7, ]
datos_alta_temperatura<- calidad_Aire_Limpio[Temperatura > 17.7, ]

length(datos_baja_temperatura$Monoxido_carbono)
## [1] 3694
length(datos_alta_temperatura$Monoxido_carbono)
## [1] 3247

Vamos a evaluar si existe una diferencia significativa en el promedio de la concentración de monóxido de carbono entre los dos grupos que hemos creado. Para hacer esto, utilizaremos un intervalo de confianza del 90%.

Inicialmente, observamos que la distribución del monóxido de carbono no sigue una distribución normal. Sin embargo, para calcular el intervalo de confianza del promedio, basta con tener mas de 25 observaciones para poder decir que los promedios se distribuyen normales. Afortunadamente, contamos con 3694 observaciones en el grupo datos_baja_temperatura y 3247 observaciones en el grupo datos_alta_temperatura, por lo que se cumple este requisito.

Dado que la varianza de la población es desconocida, utilizaremos un intervalo de confianza para saber si podemos asumir que son iguales.

# Veamos si es valido suponer varianzas iguales.
# Verifiquemos mediante un intervalo al 95%.
vartest.temp<- var.test(datos_alta_temperatura$Monoxido_carbono,
                        datos_baja_temperatura$Monoxido_carbono,
                        alternative = "two.sided",conf.level = 0.95)
vartest.temp$conf.int
## [1] 0.7410043 0.8467492
## attr(,"conf.level")
## [1] 0.95

Recordemos que con la prueba var.test lo que estamos viendo es \(\sigma_1/\sigma_2\), si ambas varianzas son iguales el intervalo de confianza de la prueba debe contener al 1.
El intervalo de confianza calculado al 95% de confianza es: \((0.741, 0.846)\), que claramente no contiene al 1, por lo tanto no podemos asumir varianzas iguales.
Finalmente podemos asumir que nuestras observaciones son independientes. Encontremos entonces el intervalo de confianza de la diferencia de promedios de la concentración de monóxido de carbono entre los dos grupos que hemos creado.

# Realizamos la prueba
prueba_promedios<- t.test(datos_baja_temperatura$Monoxido_carbono,
                          datos_alta_temperatura$Monoxido_carbono,
                          alternative="two.sided",mu=0,paired=F,var.equal=F)

prueba_promedios$conf.int
## [1] -0.09108558  0.04382863
## attr(,"conf.level")
## [1] 0.95
El intervalo de confianza al 95% para la comparación del promedio de la variable respuesta entre los dos grupos es:
\(\mu_{baja_temp} - \mu_{alta_temp} \in (-0.091, 0.043)\)


Observamos que el valor cero se encuentra dentro de este intervalo, además de ser un intervalo bastante estrecho con una longitud de \(L = 0.134\). Por lo tanto, llegamos a la conclusión de que no existen diferencias significativas en el promedio de la concentración de monóxido de carbono entre ambos grupos. Esto indica que la concentración promedio de monóxido de carbono no experimenta cambios significativos en relación a si las temperaturas registradas son superiores o inferiores a la temperatura promedio.


Ahora procederemos a desarrollar un modelo que analice la posible relación entre la concentración de monóxido de carbono y las concentraciones promedio de benceno, \(NO_x\) y \(NO_2\), así como la influencia de la temperatura y la humedad relativa en dichas concentraciones.
Antes de proceder a la formulación de un modelo, es importante explorar más a fondo las relaciones entre las distintas variables. A continuación, analizaremos estas relaciones mediante representaciones gráficas:


library(GGally)
## Warning: package 'GGally' was built under R version 4.3.2
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.3.2
## 
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
## 
##     %+%, alpha
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
ggpairs(calidad_Aire_Limpio)

Veamos una gráfica más, una representación de la matriz de correlación:


library(corrplot)
## Warning: package 'corrplot' was built under R version 4.3.2
## corrplot 0.92 loaded
corrplot(cor(calidad_Aire_Limpio))

De las gráficas presentadas podemos concluir que las variables que tienen una mayor relación lineal con la concentración de Monóxido de carboono son:

  1. Concentración de Benceno (0.930).

  2. Concentración de Óxidos de Nitrógeno (0.786).

  3. Concentración de Dióxido de Nitrógeno (0.674).

Las variables predictoras que están más relacionadas entre sí son:

  1. Óxidos de Nitrógeno - Benceno (0.718).

  2. Dióxido de Nitrógeno - Benceno (0.604).

  3. Dióxido de Nitrógeno - Óxidos de Nitrogéno (0.757).

  4. Humedad Relativa - Temperatura (0.564).

Vamos a proceder a desarrollar un primer modelo que incluya todas las variables predictoras y evaluaremos su rendimiento inicial. Sin embargo, desde un primer vistazo, podemos anticipar que el modelo no será óptimo debido a la presencia de dependencia lineal entre varias variables predictoras. Además, debemos tener en cuenta que no estamos considerando en este modelo las variables que hemos transformado previamente y que las variables predictoras no siguen una distribución normal.


Primer Modelo 

attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 6):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 7):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 9):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 10):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 11):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 12):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 13):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 15):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 21):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
# Modelo con todos los predictores
model.monoxido1<- lm(Monoxido_carbono~Benceno+Oxidos_nitrogeno+Dioxido_nitrogeno+Temperatura+Humedad_relativa, data = calidad_Aire_Limpio)  

summary(model.monoxido1)
## 
## Call:
## lm(formula = Monoxido_carbono ~ Benceno + Oxidos_nitrogeno + 
##     Dioxido_nitrogeno + Temperatura + Humedad_relativa, data = calidad_Aire_Limpio)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.2598 -0.1851  0.0062  0.1922  3.0000 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        1.836e-01  4.208e-02   4.364  1.3e-05 ***
## Benceno            1.566e-01  1.316e-03 119.015  < 2e-16 ***
## Oxidos_nitrogeno   8.093e-04  5.472e-05  14.788  < 2e-16 ***
## Dioxido_nitrogeno  2.465e-03  2.037e-04  12.099  < 2e-16 ***
## Temperatura       -1.213e-02  9.522e-04 -12.736  < 2e-16 ***
## Humedad_relativa   1.585e-03  4.351e-04   3.643 0.000271 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4529 on 6935 degrees of freedom
## Multiple R-squared:  0.9013, Adjusted R-squared:  0.9012 
## F-statistic: 1.267e+04 on 5 and 6935 DF,  p-value: < 2.2e-16

INTERPRETACION DEL PRIMER MODELO


Al parecer nuestra suposición inicial no era tan correcta, el primer modelo utilizando todas las variables predictoras no resulto tan malo, veamos algunas conclusiones:

  1. El modelo con todas las variables introducidas como predictores tiene un R alto (0.90), es capaz de explicar el 90% de la variabilidad observada en la concentración de Monóxido de Carbono.

  2. El valor-p del primer modelo es significativo (< 2.2e-16) por lo que se puede aceptar que el modelo no es por azar, al menos uno de los coeficientes parciales de regresión es distinto de 0.

  3. Todas las variables son significativas (valor-p de todos los coeficientes del modelo).

El modelo quedaría expresado como:
\(CO = 0.1836 + (0.1566)*Benceno + (0.0008)*NO_x + (0.0024)*NO_2 - (0.0123)*T + (0.0015)*Humedad\)

Validación estadística del modelo


Veamos un resumen del modelo:

library(gvlma)
valida_model <- gvlma(model.monoxido1)  
summary(valida_model)
## 
## Call:
## lm(formula = Monoxido_carbono ~ Benceno + Oxidos_nitrogeno + 
##     Dioxido_nitrogeno + Temperatura + Humedad_relativa, data = calidad_Aire_Limpio)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.2598 -0.1851  0.0062  0.1922  3.0000 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        1.836e-01  4.208e-02   4.364  1.3e-05 ***
## Benceno            1.566e-01  1.316e-03 119.015  < 2e-16 ***
## Oxidos_nitrogeno   8.093e-04  5.472e-05  14.788  < 2e-16 ***
## Dioxido_nitrogeno  2.465e-03  2.037e-04  12.099  < 2e-16 ***
## Temperatura       -1.213e-02  9.522e-04 -12.736  < 2e-16 ***
## Humedad_relativa   1.585e-03  4.351e-04   3.643 0.000271 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4529 on 6935 degrees of freedom
## Multiple R-squared:  0.9013, Adjusted R-squared:  0.9012 
## F-statistic: 1.267e+04 on 5 and 6935 DF,  p-value: < 2.2e-16
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = model.monoxido1) 
## 
##                        Value p-value                   Decision
## Global Stat        2.193e+05 0.00000 Assumptions NOT satisfied!
## Skewness           2.013e+03 0.00000 Assumptions NOT satisfied!
## Kurtosis           2.168e+05 0.00000 Assumptions NOT satisfied!
## Link Function      4.058e+00 0.04397 Assumptions NOT satisfied!
## Heteroscedasticity 4.440e+02 0.00000 Assumptions NOT satisfied!

Como ya lo habíamos comentado, todas las variables son significativas (valor-p de todos los coeficientes del modelo) y deberían quedarse en el modelo.
Los errores estándar de cada uno de los coeficientes son:

Coeficiente Error Estándar
(Intercept) 4.208e-02
Benceno 1.316e-03
Oxidos_nitrogeno 5.472e-05
Dioxido_nitrogeno 2.037e-04
Temperatura 9.522e-04
Humedad_relativa 4.351e-04

Y el error estándar residual del modelo es: 0.4529.

Diagnóstico mediante las gráficas que proporciona el programa R

# Graficas de diagnostico
plot(model.monoxido1)

Veamoos algunoos comentarios respecto a las gráficas.

  1. Gráfica residuos vs valores ajustados: Los residuos parecen estar concentrandose en la primera parte de la gráfica, mostrando alguna especie de patrón. Las observaciones 5817, 4258 y 4259 presentan diferencias mayores, por lo que quizá estas observaciones “tengan algún problema”.

  2. Normalidad residuos: Los residuos parecen no comportarse con normalidad, y vemos los mayores problemas en las “colas”.

  3. Gráfica localización-escala: Vemos que la varianza de los residuos ocomienza a crecer, formando una especie de cono. La línea roja va en aumento y no es hoorizoontal.

  4. Gráfica residuos vs leverage: El residuo 4259 está fuera de la distancia de Cook.

Lo que parecía un buen modelo por el alto valor de \(R^2\) vemos que tiene problemas muy importantes, no está cumpliendo ninguno de los supuestos para poder hacer la regresión lineal, veamos si podemos mejorarlo un poco.

Para comenzar, veamos si podemos elegir una mejor “combinación” de variable predictoras: 
step(object = model.monoxido1, direction = "both", trace = 1)
## Start:  AIC=-10989.16
## Monoxido_carbono ~ Benceno + Oxidos_nitrogeno + Dioxido_nitrogeno + 
##     Temperatura + Humedad_relativa
## 
##                     Df Sum of Sq    RSS      AIC
## <none>                           1422.6 -10989.2
## - Humedad_relativa   1      2.72 1425.3 -10977.9
## - Dioxido_nitrogeno  1     30.03 1452.6 -10846.2
## - Temperatura        1     33.28 1455.9 -10830.7
## - Oxidos_nitrogeno   1     44.86 1467.5 -10775.7
## - Benceno            1   2905.65 4328.3  -3268.1
## 
## Call:
## lm(formula = Monoxido_carbono ~ Benceno + Oxidos_nitrogeno + 
##     Dioxido_nitrogeno + Temperatura + Humedad_relativa, data = calidad_Aire_Limpio)
## 
## Coefficients:
##       (Intercept)            Benceno   Oxidos_nitrogeno  Dioxido_nitrogeno  
##         0.1836346          0.1566329          0.0008093          0.0024648  
##       Temperatura   Humedad_relativa  
##        -0.0121272          0.0015850
Como ya habíamos anticipado, según los valores p que observamos anteriormente, todas las variables resultan significativas y, por lo tanto, deben ser incluidas en el modelo. Sin embargo, con el fin de abordar el problema de la normalización, propongamos ahora la incorporación de las variables transformadas que creamos previamente. De esta manera, podremos evaluar si se produce alguna mejora en el modelo. 
# Agregamos los datos al conjunto
calidad_Aire_Limpio$transformacion_benceno <- trans.benceno
calidad_Aire_Limpio$transformacion_oxidos <- trans.oxidos
Creamos el modelo: 
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 4):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 8):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 9):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 11):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 12):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 13):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 14):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 15):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 17):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 23):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
# Modelo con todos los predictores transformados
model.monoxido2<- lm(Monoxido_carbono~transformacion_benceno+transformacion_oxidos+Dioxido_nitrogeno+Temperatura+Humedad_relativa, data = calidad_Aire_Limpio)

summary(model.monoxido2)
## 
## Call:
## lm(formula = Monoxido_carbono ~ transformacion_benceno + transformacion_oxidos + 
##     Dioxido_nitrogeno + Temperatura + Humedad_relativa, data = calidad_Aire_Limpio)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.1199 -0.3180 -0.0677  0.2559  5.4020 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            -0.5265823  0.0666031  -7.906 3.06e-15 ***
## transformacion_benceno  0.8537410  0.0103556  82.443  < 2e-16 ***
## transformacion_oxidos   0.0394882  0.0091440   4.318 1.59e-05 ***
## Dioxido_nitrogeno       0.0014394  0.0003096   4.649 3.40e-06 ***
## Temperatura            -0.0263575  0.0014257 -18.487  < 2e-16 ***
## Humedad_relativa        0.0014952  0.0005825   2.567   0.0103 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.597 on 6935 degrees of freedom
## Multiple R-squared:  0.8285, Adjusted R-squared:  0.8284 
## F-statistic:  6701 on 5 and 6935 DF,  p-value: < 2.2e-16
# Graficas de diagnostico
plot(model.monoxido2)

El modelo para nada mejoro, el valor de \(R^2\) disminuyo, la normalidad no mejoró y en general todo empeoró, lo único “bueno” es que ahora no hay observaciones fuera de la distancia de Cook.

Multicolinealidad


Ya habíamos dicho que uno de los supuestos que debe satisfacer el modelo de regresión múltiple es el de no multicolinealidad, para ver si existe usaremos: Tolerancia y Factor de Inflación de la Varianza (VIF), si:

  • VIF > 10 es causa de preocupación.

  • VIF es sustancialmente mayor que 1 entonces la regresión puede verse perjudicada.

  • Tolerancia = 1/VIF debajo de 0.1 indica un problema serio.

  • Tolerancia debajo de 0.2 indica un problema potencial.

Como el segundo modelo falló regresemos al primero. 
# Criterio VIF
library(car)
## Warning: package 'car' was built under R version 4.3.2
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.3.2
## 
## Attaching package: 'car'
## The following object is masked from 'package:psych':
## 
##     logit
## The following object is masked from 'package:dplyr':
## 
##     recode
vif(model.monoxido1)
##           Benceno  Oxidos_nitrogeno Dioxido_nitrogeno       Temperatura 
##          3.265584          4.408959          3.163125          2.399888 
##  Humedad_relativa 
##          1.946206
Observamos que efectivamente existe una correlación entre las variables, como se pudo apreciar en las gráficas anteriores que mostraban el factor de correlación. Basándonos en esta información, podemos utilizarla como criterio para eliminar algunas variables del modelo. De manera informal, hemos decidido excluir la concentración de óxidos de carbono debido a su fuerte correlación con el dióxido de carbono (0.757). Asimismo, hemos optado por eliminar la humedad relativa debido a su estrecha relación con la temperatura (-0.564).
Realicemos un nuevo modelo para ver si mejora: 
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 5):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura, transformacion_benceno,
##     transformacion_oxidos
## The following objects are masked from calidad_Aire_Limpio (pos = 7):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 11):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 12):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 14):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 15):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 16):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 17):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 18):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 20):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 26):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
# Modelo con todos los predictores transformados
model.monoxido3<- lm(Monoxido_carbono~Benceno+Dioxido_nitrogeno+Temperatura, data = calidad_Aire_Limpio)

summary(model.monoxido3)
## 
## Call:
## lm(formula = Monoxido_carbono ~ Benceno + Dioxido_nitrogeno + 
##     Temperatura, data = calidad_Aire_Limpio)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.7656 -0.1933  0.0026  0.1911  3.1530 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0.3489067  0.0213727   16.32   <2e-16 ***
## Benceno            0.1709388  0.0010206  167.49   <2e-16 ***
## Dioxido_nitrogeno  0.0034282  0.0001614   21.25   <2e-16 ***
## Temperatura       -0.0203306  0.0007031  -28.92   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4619 on 6937 degrees of freedom
## Multiple R-squared:  0.8973, Adjusted R-squared:  0.8973 
## F-statistic: 2.021e+04 on 3 and 6937 DF,  p-value: < 2.2e-16
library(gvlma)
valida_model <- gvlma(model.monoxido3)  
summary(valida_model)
## 
## Call:
## lm(formula = Monoxido_carbono ~ Benceno + Dioxido_nitrogeno + 
##     Temperatura, data = calidad_Aire_Limpio)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.7656 -0.1933  0.0026  0.1911  3.1530 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0.3489067  0.0213727   16.32   <2e-16 ***
## Benceno            0.1709388  0.0010206  167.49   <2e-16 ***
## Dioxido_nitrogeno  0.0034282  0.0001614   21.25   <2e-16 ***
## Temperatura       -0.0203306  0.0007031  -28.92   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4619 on 6937 degrees of freedom
## Multiple R-squared:  0.8973, Adjusted R-squared:  0.8973 
## F-statistic: 2.021e+04 on 3 and 6937 DF,  p-value: < 2.2e-16
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = model.monoxido3) 
## 
##                        Value p-value                   Decision
## Global Stat        2.719e+05   0.000 Assumptions NOT satisfied!
## Skewness           2.010e+03   0.000 Assumptions NOT satisfied!
## Kurtosis           2.693e+05   0.000 Assumptions NOT satisfied!
## Link Function      1.874e+00   0.171    Assumptions acceptable.
## Heteroscedasticity 5.628e+02   0.000 Assumptions NOT satisfied!
# Graficas de diagnostico
plot(model.monoxido2)

El modelo sigue sin mejorar, el primer modelo parece ser el mejor pero como ya mencionamos no cumple con los supuestos de la regresión.
Veamos ahora como úlitmo intento los valores atípicos.

Valores atípicos 

boxplot.matrix(as.matrix(calidad_Aire_Limpio),use.cols = T)
En R existe otra paqueteria que nos ayuda graficamente a encontrar estas observaciones atipicas. 
library(olsrr)
## Warning: package 'olsrr' was built under R version 4.3.2
## 
## Attaching package: 'olsrr'
## The following object is masked from 'package:MASS':
## 
##     cement
## The following object is masked from 'package:datasets':
## 
##     rivers
ols_plot_cooksd_chart(model.monoxido1)

# DFBETA mide la diferencia en la estimacion de cada parametro con y sin el punto influyente.
# Hay un valor de DFBETA para cada observacion, es decir, hay n observaciones y k parametros, por lo que hay n*k DFBETAs. 
# En general, valores grandes de DFBETAS indican observaciones que son influyentes en la estimacion de cierto parametro.

ols_plot_dfbetas(model.monoxido1)

ols_plot_resid_stud(model.monoxido1)
Como soon muchos datos la realidad es que las graficas parecen noo decir mucho, veamos entonces: 
library(car)
outlierTest(model.monoxido1)
##        rstudent unadjusted p-value Bonferroni p
## 4259 -18.844505         2.7147e-77   1.8843e-73
## 4258 -14.065135         2.5194e-44   1.7487e-40
## 5817  -7.434707         1.1740e-13   8.1486e-10
## 4918   6.648333         3.1897e-11   2.2140e-07
## 4204   6.577777         5.1231e-11   3.5559e-07
## 5418  -6.274431         3.7204e-10   2.5824e-06
## 4260  -6.111312         1.0417e-09   7.2304e-06
## 5417  -6.103445         1.0940e-09   7.5936e-06
## 5406  -5.940388         2.9811e-09   2.0692e-05
## 5407  -5.816380         6.2812e-09   4.3598e-05
Con R podemos tener un criterio para detectar outliers, descartemos entonces las observaciones de arriba y volvamos a correr el modelo 1 sin estas observaciones. 
registros_a_eliminar <- c(4259, 4258, 5817, 4918, 4204, 5418, 4260, 5417, 5406, 5407)

datos_sin_registros <- calidad_Aire_Limpio[-registros_a_eliminar, ]

model.monoxido4<- lm(Monoxido_carbono~Benceno+Oxidos_nitrogeno+Dioxido_nitrogeno+Temperatura+Humedad_relativa, data = datos_sin_registros)  
summary(model.monoxido4)
## 
## Call:
## lm(formula = Monoxido_carbono ~ Benceno + Oxidos_nitrogeno + 
##     Dioxido_nitrogeno + Temperatura + Humedad_relativa, data = datos_sin_registros)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.58037 -0.18460  0.00603  0.19348  2.56970 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        2.314e-01  3.949e-02   5.861 4.82e-09 ***
## Benceno            1.619e-01  1.250e-03 129.578  < 2e-16 ***
## Oxidos_nitrogeno   7.070e-04  5.143e-05  13.747  < 2e-16 ***
## Dioxido_nitrogeno  2.345e-03  1.910e-04  12.278  < 2e-16 ***
## Temperatura       -1.468e-02  8.972e-04 -16.366  < 2e-16 ***
## Humedad_relativa   1.270e-03  4.079e-04   3.114  0.00186 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.424 on 6925 degrees of freedom
## Multiple R-squared:  0.9131, Adjusted R-squared:  0.913 
## F-statistic: 1.455e+04 on 5 and 6925 DF,  p-value: < 2.2e-16
plot(model.monoxido4)

Después de eliminar las observaciones influyentes, hemos observado una mejora en el valor de \(R^2\) en nuestro modelo, lo cual puede parecer positivo a simple vista. Sin embargo, al examinar detenidamente las demás métricas y las últimas gráficas, es evidente que el modelo no es sólido. A pesar de que tradicionalmente utilizamos el valor de \(R^2\) como una medida para determinar la calidad del modelo, en este caso nos encontramos con problemas que incumplen los supuestos necesarios para realizar la regresión. Ninguno de estos supuestos se cumple, lo cual es motivo de preocupación y nos lleva a cuestionar si este modelo puede ser utilizado de manera apropiada.

Es posible que podríamos realizar más manipulaciones en los datos, como transformar todas las variables a una distribución normal utilizando el método de Box-Cox o explorar diferentes combinaciones de las variables predictoras según su correlación. Sin embargo, a pesar de estas opciones, concluimos que el modelo 1, a pesar de sus limitaciones y errores mencionados anteriormente, es el mejor que podemos obtener en la situación actual. Sin embargo, es importante tener en cuenta que su utilidad puede ser cuestionable y no se recomienda su uso sin tener en cuenta sus limitaciones.

Fase 4. Implementación


Con nuestro mejor modelo (modelo 1) realicemos 3 pronósticos para el valor promedio de la variable respuesta (concentración de monóxido de carbono), recordemos que nuestro modelo es:
\(CO = 0.1836 + (0.1566)*Benceno + (0.0008)*NO_x + (0.0024)*NO_2 - (0.0123)*T + (0.0015)*Humedad\) 
attach(calidad_Aire_Limpio)
## The following objects are masked from calidad_Aire_Limpio (pos = 4):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura, transformacion_benceno,
##     transformacion_oxidos
## The following objects are masked from calidad_Aire_Limpio (pos = 7):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura, transformacion_benceno,
##     transformacion_oxidos
## The following objects are masked from calidad_Aire_Limpio (pos = 9):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 13):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 14):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 16):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 17):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 18):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 19):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 20):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 22):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
## The following objects are masked from calidad_Aire_Limpio (pos = 28):
## 
##     Benceno, Dioxido_nitrogeno, Humedad_relativa, Monoxido_carbono,
##     Oxidos_nitrogeno, Temperatura
nuevos_datos <- data.frame(
  var1 = c(max(Benceno)+5,
           max(Oxidos_nitrogeno)+20,
           max(Dioxido_nitrogeno)+20,
           max(Temperatura)+3,
           max(Humedad_relativa)+2),  # Valores para var1
  var2 = c(max(Benceno)+10,
           max(Oxidos_nitrogeno)+40,
           max(Dioxido_nitrogeno)+40,
           max(Temperatura)+5,
           max(Humedad_relativa)+4),  # Valores para var2
  var3 = c(max(Benceno)+15,
           max(Oxidos_nitrogeno)+50,
           max(Dioxido_nitrogeno)+50,
           max(Temperatura)+7,
           max(Humedad_relativa)+6),  # Valores para var3
  stringsAsFactors = FALSE
)

predicciones <- predict(model.monoxido1, nuevos_datos, interval="prediction")
## Warning: 'newdata' had 5 rows but variables found have 6941 rows
head(predicciones)
##        fit       lwr      upr
## 1 2.370100 1.4820364 3.258164
## 2 1.879968 0.9918811 2.768055
## 3 1.921226 1.0331577 2.809294
## 4 2.030767 1.1426972 2.918838
## 5 1.555741 0.6676693 2.443813
## 6 1.193149 0.3050997 2.081198